      subroutine bec8g2(coorr,coefr,
     & prmt,estif,emass,edamp,eload,num)
c .... coorr ---- nodal coordinate value
c .... coefr ---- nodal coef value
      implicit real*8 (a-h,o-z)
      dimension estif(64,64),elump(64),emass(64),
     & eload(64)
      dimension prmt(*),coef(3),coefr(8,3),coorr(3,8),coor(3)
      common /rbec8g2/rsxx(8,32),rsyy(8,32),rszz(8,32),
     & rsyz(8,32),rsxz(8,32),rsxy(8,32),rsq(8,32),
     & rdid(8,32),
     & csxx(8,4),csyy(8,4),cszz(8,4),csyz(8,4),
     & csxz(8,4),csxy(8,4),csq(8,4),cdid(8,4)
c .... store shape functions and their partial derivatives
c .... for all integral points
      common /vbec8g2/rctr(3,3),crtr(3,3),coefd(3,9),coefc(3,9)
      common /dbec8g2/ refc(3,8),gaus(8),
     & nnode,ngaus,ndisp,nrefc,ncoor,nvar,
     & nvard(8),kdord(8),kvord(64,8)
c .... nnode ---- the number of nodes
c .... nrefc ---- the number of numerical integral points
c .... ndisp ---- the number of unknown functions
c .... nrefc ---- the number of reference coordinates
c .... nvar ---- the number of unknown varibles var
c .... refc ---- reference coordinates at integral points
c .... gaus ---- weight number at integral points
c .... nvard ---- the number of var for each unknown
c .... kdord ---- the highest differential order for each unknown
c .... kvord ---- var number at integral points for each unknown
      dimension dev(3),dep(3),ds(6),dspv(3),dspp(3),
     & dpv(3),dpp(3),e(3,3),d(3,3)
      dimension ddf(6),p(4)
      external amohrc
      common /gpstr/ gpstr(7,1000000),str(7)
      nstr = 7
      pe=prmt(1)
      pv=prmt(2)
      fx=prmt(3)
      fy=prmt(4)
      fz=prmt(5)
      p(1)=prmt(6)
      p(2)=prmt(7)
      p(4)=prmt(8)
      qnn=prmt(9)
      rou=prmt(10)
      wrou=prmt(11)
      grndwt_lv=prmt(12)
      alpha=prmt(13)
      edjt=prmt(14)
      time=prmt(15)
      dt=prmt(16)
      imate=prmt(17)+0.5
      ielem=prmt(18)+0.5
      nelem=prmt(19)+0.5
      it=prmt(20)+0.5
      nmate=prmt(21)+0.5
      itime=prmt(22)+0.5
      ityp=prmt(23)+0.5
      kq=int(qnn+0.5)
      fact = pe/(1.+pv)/(1.-pv*2)
      shear = (0.5-pv)
      p(3)=pv
      shear=shear*fact
      if (num.eq.1) call bec8g2i
c .... initialize the basic data
      do 10 i=1,nvar
      emass(i)=0.0
      eload(i)=0.0
      do 10 j=1,nvar
      estif(i,j)=0.0
10    continue
      do 999 igaus=1,ngaus
      call bec8g2t(nnode,nrefc,ncoor,refc(1,igaus),coor,coorr,
     & rctr,crtr,det,coefr)
c .... coordinate transfer from reference to original system
c .... rctr ---- Jacobi's matrix
c .... crtr ---- inverse matrix of Jacobi's matrix
      x=coor(1)
      y=coor(2)
      z=coor(3)
      rx=refc(1,igaus)
      ry=refc(2,igaus)
      rz=refc(3,igaus)
      call ebec8g2(refc(1,igaus),coef,coorr,coefr,coefd)
c .... compute coef functions and their partial derivatives
      isxx=(igaus-1)*4+1
      isyy=(igaus-1)*4+1
      iszz=(igaus-1)*4+1
      isyz=(igaus-1)*4+1
      isxz=(igaus-1)*4+1
      isxy=(igaus-1)*4+1
      isq=(igaus-1)*4+1
      idid=(igaus-1)*4+1
      if (num.gt.1) goto 2
c .... the following is the shape function caculation
      call bec8g21(refc(1,igaus),rsxx(1,isxx),rctr,crtr)
      call bec8g22(refc(1,igaus),rsyy(1,isyy),rctr,crtr)
      call bec8g23(refc(1,igaus),rszz(1,iszz),rctr,crtr)
      call bec8g24(refc(1,igaus),rsyz(1,isyz),rctr,crtr)
      call bec8g25(refc(1,igaus),rsxz(1,isxz),rctr,crtr)
      call bec8g26(refc(1,igaus),rsxy(1,isxy),rctr,crtr)
      call bec8g27(refc(1,igaus),rsq(1,isq),rctr,crtr)
      call bec8g28(refc(1,igaus),rdid(1,idid),rctr,crtr)
2     continue
c .... the following is the shape function transformation
c .... from reference coordinates to original coordinates
      call shapn(nrefc,ncoor,8,rsxx(1,isxx),csxx,crtr,1,4,4)
      call shapn(nrefc,ncoor,8,rsyy(1,isyy),csyy,crtr,1,4,4)
      call shapn(nrefc,ncoor,8,rszz(1,iszz),cszz,crtr,1,4,4)
      call shapn(nrefc,ncoor,8,rsyz(1,isyz),csyz,crtr,1,4,4)
      call shapn(nrefc,ncoor,8,rsxz(1,isxz),csxz,crtr,1,4,4)
      call shapn(nrefc,ncoor,8,rsxy(1,isxy),csxy,crtr,1,4,4)
      call shapn(nrefc,ncoor,8,rsq(1,isq),csq,crtr,1,4,4)
      call shapn(nrefc,ncoor,8,rdid(1,idid),cdid,crtr,1,4,4)
c .... the coef function transformation
c .... from reference coordinates to original coordinates
      call shapc(nrefc,ncoor,3,coefd,coefc,crtr,2,9,9)
      un=coef(1)
      vn=coef(2)
      wn=coef(3)
      weigh=det*gaus(igaus)
      d(1,1) = +(1.-pv)*fact
      d(1,2) = +pv*fact
      d(1,3) = +pv*fact
      d(2,1) = +pv*fact
      d(2,2) = +(1.-pv)*fact
      d(2,3) = +pv*fact
      d(3,1) = +pv*fact
      d(3,2) = +pv*fact
      d(3,3) = +(1.-pv)*fact
      vol = 1.0
      if ((time.gt.1.5*dt)
     c .and. (num.eq.1) .and. (igaus.eq.1)) then
      open(17,file='filestr',form='unformatted',status='old')
      read(17) ((gpstr(i,j),i=1,nstr),j=1,nelem*ngaus)
      close(17)
      endif
      if (time.gt.1.5*dt) then
      do i=1,nstr
      str(i)=gpstr(i,(num-1)*ngaus+igaus)
      enddo
      else
      do i=1,nstr
      str(i)=0.0d0
      enddo
      endif
      e(1,1)=coefc(1,1)
      e(1,2)=coefc(1,2)
      e(1,3)=coefc(1,3)
      e(2,1)=coefc(2,1)
      e(2,2)=coefc(2,2)
      e(2,3)=coefc(2,3)
      e(3,1)=coefc(3,1)
      e(3,2)=coefc(3,2)
      e(3,3)=coefc(3,3)
      dev(1) = +e(1,1)
      dev(2) = +e(2,2)
      dev(3) = +e(3,3)
      dep(1) = e(2,3)+e(3,2)
      dep(2) = e(1,3)+e(3,1)
      dep(3) = e(1,2)+e(2,1)
      ds(1) = +d(1,1)*dev(1)+d(1,2)*dev(2)+d(1,3)*dev(3)
      ds(2) = +d(2,1)*dev(1)+d(2,2)*dev(2)+d(2,3)*dev(3)
      ds(3) = +d(3,1)*dev(1)+d(3,2)*dev(2)+d(3,3)*dev(3)
      ds(4) = shear*dep(1)
      ds(5) = shear*dep(2)
      ds(6) = shear*dep(3)
      ialpha=0
      prag=amohrc(p,str)
      str(1) = str(1)+ds(1)
      str(2) = str(2)+ds(2)
      str(3) = str(3)+ds(3)
      str(4) = str(4)+ds(4)
      str(5) = str(5)+ds(5)
      str(6) = str(6)+ds(6)
      call getyield(p,str,d,shear,a,ddf,kq,amohrc,ialpha)
      prag1=amohrc(p,str)
      dln=+ddf(1)*dev(1)+ddf(4)*dep(1)+ddf(2)*dev(2)
     * +ddf(5)*dep(2)+ddf(3)*dev(3)+ddf(6)*dep(3)
      dln=dln+prag
      if (dln.gt.1.0d-6) then
c	    if(prag1.gt.-p(2)*0.0001) then
       b=1.0/a
       eid=1.0d0
      ds(1) = +b*ddf(1)*ddf(1)*dev(1)+b*ddf(1)*ddf(2)
     * *dev(2)+b*ddf(1)*ddf(3)*dev(3)
      ds(2) = +b*ddf(2)*ddf(1)*dev(1)+b*ddf(2)*ddf(2)
     * *dev(2)+b*ddf(2)*ddf(3)*dev(3)
      ds(3) = +b*ddf(3)*ddf(1)*dev(1)+b*ddf(3)*ddf(2)
     * *dev(2)+b*ddf(3)*ddf(3)*dev(3)
      ds(4) = +b*ddf(4)*ddf(4)*dep(1)
      ds(5) = +b*ddf(5)*ddf(5)*dep(2)
      ds(6) = +b*ddf(6)*ddf(6)*dep(3)
      str(1) = str(1)-ds(1)-ddf(1)*prag*b
      str(2) = str(2)-ds(2)-ddf(2)*prag*b
      str(3) = str(3)-ds(3)-ddf(3)*prag*b
      str(4) = str(4)-ds(4)-ddf(4)*prag*b
      str(5) = str(5)-ds(5)-ddf(5)*prag*b
      str(6) = str(6)-ds(6)-ddf(6)*prag*b
      else
       b=0.0d0
       eid=0.0d0
      end if
      call getyield(p,str,d,shear,a,ddf,kq,amohrc,ialpha)
      cv=+ddf(1)*dev(1)+ddf(2)*dev(2)+ddf(3)*dev(3)
      cp=+ddf(4)*dep(1)+ddf(5)*dep(2)+ddf(6)*dep(3)
      dspv(1) = ddf(1)*(cv+cp)*b
      dspv(2) = ddf(2)*(cv+cp)*b
      dspv(3) = ddf(3)*(cv+cp)*b
      dspp(1) = ddf(4)*(cv+cp)*b
      dspp(2) = ddf(5)*(cv+cp)*b
      dspp(3) = ddf(6)*(cv+cp)*b
      detd = 0.0
      call invm(3,d,detd)
      dpv(1) = +d(1,1)*dspv(1)+d(1,2)*dspv(2)+d(1,3)
     * *dspv(3)
      dpv(2) = +d(2,1)*dspv(1)+d(2,2)*dspv(2)+d(2,3)
     * *dspv(3)
      dpv(3) = +d(3,1)*dspv(1)+d(3,2)*dspv(2)+d(3,3)
     * *dspv(3)
      dpp(1) = dspp(1)/shear
      dpp(2) = dspp(2)/shear
      dpp(3) = dspp(3)/shear
      if (kq.eq.2) then
      epp=+dpv(1)+dpv(2)+dpv(3)
      str(7)=epp+str(7)
      elseif (kq.eq.3) then
      eyp=+dpv(1)**2+dpp(1)**2+dpv(2)**2+dpp(2)**2+dpv(3)
     * **2+dpp(3)**2
      str(7)=dsqrt(eyp)+str(7)
      else
      wpp=+dpv(1)*str(1)+dpp(1)*str(4)+dpv(2)*str(2)
     * +dpp(2)*str(5)+dpv(3)*str(3)+dpp(3)*str(6)
      str(7)=wpp+str(7)
      end if
      do i=1,nstr
      gpstr(i,(num-1)*ngaus+igaus)=str(i)
      enddo
      if ((num.eq.nelem) .and. (igaus.eq.ngaus)) then
      open(17,file='filestr',form='unformatted',status='unknown')
      write(17) ((gpstr(i,j),i=1,nstr),j=1,nelem*ngaus)
      close(17)
      open(21,file='filestr0',form='formatted',status='unknown')
      do j=1,nelem
        do ig=1,ngaus
          k=(j-1)*ngaus+ig
          write(21,'(2i8,8e15.7)') j,ig,(gpstr(i,k),i=1,nstr)
        enddo
      enddo
      endif
c .... the following is the stiffness computation
      do 202 i=1,8
      iv=kvord(i,1)
      do 201 j=1,8
      jv=kvord(j,1)
      stif=+csxx(i,1)*csxx(j,1)*0.0
      estif(iv,jv)=estif(iv,jv)+stif*weigh
201    continue
202    continue
c .... the following is the mass matrix computation
      stif=vol
      elump(1)=stif*weigh
      stif=vol
      elump(9)=stif*weigh
      stif=vol
      elump(17)=stif*weigh
      stif=vol
      elump(25)=stif*weigh
      stif=vol
      elump(33)=stif*weigh
      stif=vol
      elump(41)=stif*weigh
      stif=vol
      elump(49)=stif*weigh
      stif=vol
      elump(57)=stif*weigh
      stif=vol
      elump(2)=stif*weigh
      stif=vol
      elump(10)=stif*weigh
      stif=vol
      elump(18)=stif*weigh
      stif=vol
      elump(26)=stif*weigh
      stif=vol
      elump(34)=stif*weigh
      stif=vol
      elump(42)=stif*weigh
      stif=vol
      elump(50)=stif*weigh
      stif=vol
      elump(58)=stif*weigh
      stif=vol
      elump(3)=stif*weigh
      stif=vol
      elump(11)=stif*weigh
      stif=vol
      elump(19)=stif*weigh
      stif=vol
      elump(27)=stif*weigh
      stif=vol
      elump(35)=stif*weigh
      stif=vol
      elump(43)=stif*weigh
      stif=vol
      elump(51)=stif*weigh
      stif=vol
      elump(59)=stif*weigh
      stif=vol
      elump(4)=stif*weigh
      stif=vol
      elump(12)=stif*weigh
      stif=vol
      elump(20)=stif*weigh
      stif=vol
      elump(28)=stif*weigh
      stif=vol
      elump(36)=stif*weigh
      stif=vol
      elump(44)=stif*weigh
      stif=vol
      elump(52)=stif*weigh
      stif=vol
      elump(60)=stif*weigh
      stif=vol
      elump(5)=stif*weigh
      stif=vol
      elump(13)=stif*weigh
      stif=vol
      elump(21)=stif*weigh
      stif=vol
      elump(29)=stif*weigh
      stif=vol
      elump(37)=stif*weigh
      stif=vol
      elump(45)=stif*weigh
      stif=vol
      elump(53)=stif*weigh
      stif=vol
      elump(61)=stif*weigh
      stif=vol
      elump(6)=stif*weigh
      stif=vol
      elump(14)=stif*weigh
      stif=vol
      elump(22)=stif*weigh
      stif=vol
      elump(30)=stif*weigh
      stif=vol
      elump(38)=stif*weigh
      stif=vol
      elump(46)=stif*weigh
      stif=vol
      elump(54)=stif*weigh
      stif=vol
      elump(62)=stif*weigh
      stif=vol
      elump(7)=stif*weigh
      stif=vol
      elump(15)=stif*weigh
      stif=vol
      elump(23)=stif*weigh
      stif=vol
      elump(31)=stif*weigh
      stif=vol
      elump(39)=stif*weigh
      stif=vol
      elump(47)=stif*weigh
      stif=vol
      elump(55)=stif*weigh
      stif=vol
      elump(63)=stif*weigh
      stif=vol
      elump(8)=stif*weigh
      stif=vol
      elump(16)=stif*weigh
      stif=vol
      elump(24)=stif*weigh
      stif=vol
      elump(32)=stif*weigh
      stif=vol
      elump(40)=stif*weigh
      stif=vol
      elump(48)=stif*weigh
      stif=vol
      elump(56)=stif*weigh
      stif=vol
      elump(64)=stif*weigh
      do 301 i=1,nvard(1)
      iv = kvord(i,1)
      emass(iv)=emass(iv)+elump(iv)*csxx(i,1)
301   continue
      do 302 i=1,nvard(2)
      iv = kvord(i,2)
      emass(iv)=emass(iv)+elump(iv)*csyy(i,1)
302   continue
      do 303 i=1,nvard(3)
      iv = kvord(i,3)
      emass(iv)=emass(iv)+elump(iv)*cszz(i,1)
303   continue
      do 304 i=1,nvard(4)
      iv = kvord(i,4)
      emass(iv)=emass(iv)+elump(iv)*csyz(i,1)
304   continue
      do 305 i=1,nvard(5)
      iv = kvord(i,5)
      emass(iv)=emass(iv)+elump(iv)*csxz(i,1)
305   continue
      do 306 i=1,nvard(6)
      iv = kvord(i,6)
      emass(iv)=emass(iv)+elump(iv)*csxy(i,1)
306   continue
      do 307 i=1,nvard(7)
      iv = kvord(i,7)
      emass(iv)=emass(iv)+elump(iv)*csq(i,1)
307   continue
      do 308 i=1,nvard(8)
      iv = kvord(i,8)
      emass(iv)=emass(iv)+elump(iv)*cdid(i,1)
308   continue
c .... the following is the load vector computation
      do 501 i=1,8
      iv=kvord(i,1)
      stif=+csxx(i,1)*str(1)*vol
      eload(iv)=eload(iv)+stif*weigh
501   continue
      do 502 i=1,8
      iv=kvord(i,2)
      stif=+csyy(i,1)*str(2)*vol
      eload(iv)=eload(iv)+stif*weigh
502   continue
      do 503 i=1,8
      iv=kvord(i,3)
      stif=+cszz(i,1)*str(3)*vol
      eload(iv)=eload(iv)+stif*weigh
503   continue
      do 504 i=1,8
      iv=kvord(i,4)
      stif=+csyz(i,1)*str(4)*vol
      eload(iv)=eload(iv)+stif*weigh
504   continue
      do 505 i=1,8
      iv=kvord(i,5)
      stif=+csxz(i,1)*str(5)*vol
      eload(iv)=eload(iv)+stif*weigh
505   continue
      do 506 i=1,8
      iv=kvord(i,6)
      stif=+csxy(i,1)*str(6)*vol
      eload(iv)=eload(iv)+stif*weigh
506   continue
      do 507 i=1,8
      iv=kvord(i,7)
      stif=+csq(i,1)*str(7)*vol
      eload(iv)=eload(iv)+stif*weigh
507   continue
      do 508 i=1,8
      iv=kvord(i,8)
      stif=+cdid(i,1)*eid*vol
      eload(iv)=eload(iv)+stif*weigh
508   continue
999   continue
998   continue
      return
      end

      subroutine bec8g2i
      implicit real*8 (a-h,o-z)
      common /dbec8g2/ refc(3,8),gaus(8),
     & nnode,ngaus,ndisp,nrefc,ncoor,nvar,
     & nvard(8),kdord(8),kvord(64,8)
c .... initial data
c .... refc ---- reference coordinates at integral points
c .... gaus ---- weight number at integral points
c .... nvard ---- the number of var for each unknown
c .... kdord ---- the highest differential order for each unknown
c .... kvord ---- var number at integral points for each unknown
      ngaus=  8
      ndisp=  8
      nrefc=  3
      ncoor=  3
      nvar = 64
      nnode=  8
      kdord(1)=1
      nvard(1)=8
      kvord(1,1)=1
      kvord(2,1)=9
      kvord(3,1)=17
      kvord(4,1)=25
      kvord(5,1)=33
      kvord(6,1)=41
      kvord(7,1)=49
      kvord(8,1)=57
      kdord(2)=1
      nvard(2)=8
      kvord(1,2)=2
      kvord(2,2)=10
      kvord(3,2)=18
      kvord(4,2)=26
      kvord(5,2)=34
      kvord(6,2)=42
      kvord(7,2)=50
      kvord(8,2)=58
      kdord(3)=1
      nvard(3)=8
      kvord(1,3)=3
      kvord(2,3)=11
      kvord(3,3)=19
      kvord(4,3)=27
      kvord(5,3)=35
      kvord(6,3)=43
      kvord(7,3)=51
      kvord(8,3)=59
      kdord(4)=1
      nvard(4)=8
      kvord(1,4)=4
      kvord(2,4)=12
      kvord(3,4)=20
      kvord(4,4)=28
      kvord(5,4)=36
      kvord(6,4)=44
      kvord(7,4)=52
      kvord(8,4)=60
      kdord(5)=1
      nvard(5)=8
      kvord(1,5)=5
      kvord(2,5)=13
      kvord(3,5)=21
      kvord(4,5)=29
      kvord(5,5)=37
      kvord(6,5)=45
      kvord(7,5)=53
      kvord(8,5)=61
      kdord(6)=1
      nvard(6)=8
      kvord(1,6)=6
      kvord(2,6)=14
      kvord(3,6)=22
      kvord(4,6)=30
      kvord(5,6)=38
      kvord(6,6)=46
      kvord(7,6)=54
      kvord(8,6)=62
      kdord(7)=1
      nvard(7)=8
      kvord(1,7)=7
      kvord(2,7)=15
      kvord(3,7)=23
      kvord(4,7)=31
      kvord(5,7)=39
      kvord(6,7)=47
      kvord(7,7)=55
      kvord(8,7)=63
      kdord(8)=1
      nvard(8)=8
      kvord(1,8)=8
      kvord(2,8)=16
      kvord(3,8)=24
      kvord(4,8)=32
      kvord(5,8)=40
      kvord(6,8)=48
      kvord(7,8)=56
      kvord(8,8)=64
      refc(1,1)=5.773502692e-001
      refc(2,1)=5.773502692e-001
      refc(3,1)=5.773502692e-001
      gaus(1)=1.000000000e+000
      refc(1,2)=5.773502692e-001
      refc(2,2)=5.773502692e-001
      refc(3,2)=-5.773502692e-001
      gaus(2)=1.000000000e+000
      refc(1,3)=5.773502692e-001
      refc(2,3)=-5.773502692e-001
      refc(3,3)=5.773502692e-001
      gaus(3)=1.000000000e+000
      refc(1,4)=5.773502692e-001
      refc(2,4)=-5.773502692e-001
      refc(3,4)=-5.773502692e-001
      gaus(4)=1.000000000e+000
      refc(1,5)=-5.773502692e-001
      refc(2,5)=5.773502692e-001
      refc(3,5)=5.773502692e-001
      gaus(5)=1.000000000e+000
      refc(1,6)=-5.773502692e-001
      refc(2,6)=5.773502692e-001
      refc(3,6)=-5.773502692e-001
      gaus(6)=1.000000000e+000
      refc(1,7)=-5.773502692e-001
      refc(2,7)=-5.773502692e-001
      refc(3,7)=5.773502692e-001
      gaus(7)=1.000000000e+000
      refc(1,8)=-5.773502692e-001
      refc(2,8)=-5.773502692e-001
      refc(3,8)=-5.773502692e-001
      gaus(8)=1.000000000e+000
      end


      subroutine bec8g2t(nnode,nrefc,ncoor,refc,coor,coorr,
     & rc,cr,det,coefr)
      implicit real*8 (a-h,o-z)
      dimension refc(nrefc),rc(ncoor,nrefc),cr(nrefc,ncoor),a(5,10),
     & coorr(ncoor,nnode),coor(ncoor),coefr(nnode,*)
      call tbec8g2(refc,coor,coorr,coefr,rc)
      n=nrefc
      m=n*2
      det = 1.0
      do 10 i=1,n
      do 10 j=1,n
      if (i.le.ncoor) a(i,j) = rc(i,j)
      if (i.gt.ncoor) a(i,j)=1.0
      a(i,n+j)=0.0
      if (i.eq.j) a(i,n+i) = 1.0
10    continue
c     write(*,*) 'a ='
c     do 21 i=1,n
c21   write(*,8) (a(i,j),j=1,m)
      do 400 i=1,n
      amax = 0.0
      l = 0
      do 50 j=i,n
      c = a(j,i)
      if (c.lt.0.0) c = -c
      if (c.le.amax) goto 50
      amax = c
      l = j
50    continue
      do 60 k=1,m
      c = a(l,k)
      a(l,k) = a(i,k)
      a(i,k) = c
60    continue
      c = a(i,i)
      det = c*det
      do 100 k=i+1,m
100   a(i,k) = a(i,k)/c
      do 300 j=1,n
      if (i.eq.j) goto 300
      do 200 k=i+1,m
200   a(j,k) = a(j,k)-a(i,k)*a(j,i)
c     write(*,*) 'i =',i,'  j =',j,'  a ='
c     do 11 ii=1,n
c11   write(*,8) (a(ii,jj),jj=1,m)
300   continue
400   continue
      do 500 i=1,nrefc
      do 500 j=1,ncoor
500   cr(i,j) = a(i,n+j)
c     write(*,*) 'a ='
c     do 22 i=1,n
c22   write(*,8) (a(i,j),j=1,m)
c     write(*,*) 'rc ='
c     do 24 i=1,ncoor
c24   write(*,8) (rc(i,j),j=1,nrefc)
c     write(*,*) 'cr ='
c     do 23 i=1,nrefc
c23   write(*,8) (cr(i,j),j=1,ncoor)
c     write(*,*) 'det =',det
      if (det.lt.0.0) det=-det
c     write(*,*) 'det =',det
8     format(1x,6f12.3)
      end

      subroutine bec8g21(refc,shpr,rctr,crtr)
c .... compute shape functions and their partial derivatives
c .... shapr ---- store shape functions and their partial derivatives
      implicit real*8 (a-h,o-z)
      dimension refc(3),shpr(8,4),rctr(3,3),crtr(3,3)
      external fbec8g21
      rx=refc(1)
      ry=refc(2)
      rz=refc(3)
      call dshap(fbec8g21,refc,shpr,3,8,1)
c .... shape function and their derivatives computation
c .... compute partial derivatives by centered difference
c .... which is in the file ccshap.for of FEPG library
      return
      end

      real*8 function fbec8g21(refc,n)
c .... shape function caculation
      implicit real*8 (a-h,o-z)
      common /ccbec8g2/ xa(8),ya(8),za(8),una(8),
     &vna(8),wna(8)
      common /vbec8g2/ rctr(3,3),crtr(3,3),coefd(3,9),coefc(3,9)
      dimension refc(3)
      common /coord/ coor(3),coora(27,3)
      x=coor(1)
      y=coor(2)
      z=coor(3)
      rx=refc(1)
      ry=refc(2)
      rz=refc(3)
      goto (1,2,3,4,5,6,7,8) n
1     fbec8g21=+(+1.-rx)/2.*(+1.-ry)/2.*(+1.-rz)/2. 
      goto 1000
2     fbec8g21=+(+1.+rx)/2.*(+1.-ry)/2.*(+1.-rz)/2. 
      goto 1000
3     fbec8g21=+(+1.+rx)/2.*(+1.+ry)/2.*(+1.-rz)/2. 
      goto 1000
4     fbec8g21=+(+1.-rx)/2.*(+1.+ry)/2.*(+1.-rz)/2. 
      goto 1000
5     fbec8g21=+(+1.-rx)/2.*(+1.-ry)/2.*(+1.+rz)/2. 
      goto 1000
6     fbec8g21=+(+1.+rx)/2.*(+1.-ry)/2.*(+1.+rz)/2. 
      goto 1000
7     fbec8g21=+(+1.+rx)/2.*(+1.+ry)/2.*(+1.+rz)/2. 
      goto 1000
8     fbec8g21=+(+1.-rx)/2.*(+1.+ry)/2.*(+1.+rz)/2. 
      goto 1000
1000  return
      end

      subroutine bec8g22(refc,shpr,rctr,crtr)
c .... compute shape functions and their partial derivatives
c .... shapr ---- store shape functions and their partial derivatives
      implicit real*8 (a-h,o-z)
      dimension refc(3),shpr(8,4),rctr(3,3),crtr(3,3)
      external fbec8g22
      rx=refc(1)
      ry=refc(2)
      rz=refc(3)
      call dshap(fbec8g22,refc,shpr,3,8,1)
c .... shape function and their derivatives computation
c .... compute partial derivatives by centered difference
c .... which is in the file ccshap.for of FEPG library
      return
      end

      real*8 function fbec8g22(refc,n)
c .... shape function caculation
      implicit real*8 (a-h,o-z)
      common /ccbec8g2/ xa(8),ya(8),za(8),una(8),
     &vna(8),wna(8)
      common /vbec8g2/ rctr(3,3),crtr(3,3),coefd(3,9),coefc(3,9)
      dimension refc(3)
      common /coord/ coor(3),coora(27,3)
      x=coor(1)
      y=coor(2)
      z=coor(3)
      rx=refc(1)
      ry=refc(2)
      rz=refc(3)
      goto (1,2,3,4,5,6,7,8) n
1     fbec8g22=+(+1.-rx)/2.*(+1.-ry)/2.*(+1.-rz)/2. 
      goto 1000
2     fbec8g22=+(+1.+rx)/2.*(+1.-ry)/2.*(+1.-rz)/2. 
      goto 1000
3     fbec8g22=+(+1.+rx)/2.*(+1.+ry)/2.*(+1.-rz)/2. 
      goto 1000
4     fbec8g22=+(+1.-rx)/2.*(+1.+ry)/2.*(+1.-rz)/2. 
      goto 1000
5     fbec8g22=+(+1.-rx)/2.*(+1.-ry)/2.*(+1.+rz)/2. 
      goto 1000
6     fbec8g22=+(+1.+rx)/2.*(+1.-ry)/2.*(+1.+rz)/2. 
      goto 1000
7     fbec8g22=+(+1.+rx)/2.*(+1.+ry)/2.*(+1.+rz)/2. 
      goto 1000
8     fbec8g22=+(+1.-rx)/2.*(+1.+ry)/2.*(+1.+rz)/2. 
      goto 1000
1000  return
      end

      subroutine bec8g23(refc,shpr,rctr,crtr)
c .... compute shape functions and their partial derivatives
c .... shapr ---- store shape functions and their partial derivatives
      implicit real*8 (a-h,o-z)
      dimension refc(3),shpr(8,4),rctr(3,3),crtr(3,3)
      external fbec8g23
      rx=refc(1)
      ry=refc(2)
      rz=refc(3)
      call dshap(fbec8g23,refc,shpr,3,8,1)
c .... shape function and their derivatives computation
c .... compute partial derivatives by centered difference
c .... which is in the file ccshap.for of FEPG library
      return
      end

      real*8 function fbec8g23(refc,n)
c .... shape function caculation
      implicit real*8 (a-h,o-z)
      common /ccbec8g2/ xa(8),ya(8),za(8),una(8),
     &vna(8),wna(8)
      common /vbec8g2/ rctr(3,3),crtr(3,3),coefd(3,9),coefc(3,9)
      dimension refc(3)
      common /coord/ coor(3),coora(27,3)
      x=coor(1)
      y=coor(2)
      z=coor(3)
      rx=refc(1)
      ry=refc(2)
      rz=refc(3)
      goto (1,2,3,4,5,6,7,8) n
1     fbec8g23=+(+1.-rx)/2.*(+1.-ry)/2.*(+1.-rz)/2. 
      goto 1000
2     fbec8g23=+(+1.+rx)/2.*(+1.-ry)/2.*(+1.-rz)/2. 
      goto 1000
3     fbec8g23=+(+1.+rx)/2.*(+1.+ry)/2.*(+1.-rz)/2. 
      goto 1000
4     fbec8g23=+(+1.-rx)/2.*(+1.+ry)/2.*(+1.-rz)/2. 
      goto 1000
5     fbec8g23=+(+1.-rx)/2.*(+1.-ry)/2.*(+1.+rz)/2. 
      goto 1000
6     fbec8g23=+(+1.+rx)/2.*(+1.-ry)/2.*(+1.+rz)/2. 
      goto 1000
7     fbec8g23=+(+1.+rx)/2.*(+1.+ry)/2.*(+1.+rz)/2. 
      goto 1000
8     fbec8g23=+(+1.-rx)/2.*(+1.+ry)/2.*(+1.+rz)/2. 
      goto 1000
1000  return
      end

      subroutine bec8g24(refc,shpr,rctr,crtr)
c .... compute shape functions and their partial derivatives
c .... shapr ---- store shape functions and their partial derivatives
      implicit real*8 (a-h,o-z)
      dimension refc(3),shpr(8,4),rctr(3,3),crtr(3,3)
      external fbec8g24
      rx=refc(1)
      ry=refc(2)
      rz=refc(3)
      call dshap(fbec8g24,refc,shpr,3,8,1)
c .... shape function and their derivatives computation
c .... compute partial derivatives by centered difference
c .... which is in the file ccshap.for of FEPG library
      return
      end

      real*8 function fbec8g24(refc,n)
c .... shape function caculation
      implicit real*8 (a-h,o-z)
      common /ccbec8g2/ xa(8),ya(8),za(8),una(8),
     &vna(8),wna(8)
      common /vbec8g2/ rctr(3,3),crtr(3,3),coefd(3,9),coefc(3,9)
      dimension refc(3)
      common /coord/ coor(3),coora(27,3)
      x=coor(1)
      y=coor(2)
      z=coor(3)
      rx=refc(1)
      ry=refc(2)
      rz=refc(3)
      goto (1,2,3,4,5,6,7,8) n
1     fbec8g24=+(+1.-rx)/2.*(+1.-ry)/2.*(+1.-rz)/2. 
      goto 1000
2     fbec8g24=+(+1.+rx)/2.*(+1.-ry)/2.*(+1.-rz)/2. 
      goto 1000
3     fbec8g24=+(+1.+rx)/2.*(+1.+ry)/2.*(+1.-rz)/2. 
      goto 1000
4     fbec8g24=+(+1.-rx)/2.*(+1.+ry)/2.*(+1.-rz)/2. 
      goto 1000
5     fbec8g24=+(+1.-rx)/2.*(+1.-ry)/2.*(+1.+rz)/2. 
      goto 1000
6     fbec8g24=+(+1.+rx)/2.*(+1.-ry)/2.*(+1.+rz)/2. 
      goto 1000
7     fbec8g24=+(+1.+rx)/2.*(+1.+ry)/2.*(+1.+rz)/2. 
      goto 1000
8     fbec8g24=+(+1.-rx)/2.*(+1.+ry)/2.*(+1.+rz)/2. 
      goto 1000
1000  return
      end

      subroutine bec8g25(refc,shpr,rctr,crtr)
c .... compute shape functions and their partial derivatives
c .... shapr ---- store shape functions and their partial derivatives
      implicit real*8 (a-h,o-z)
      dimension refc(3),shpr(8,4),rctr(3,3),crtr(3,3)
      external fbec8g25
      rx=refc(1)
      ry=refc(2)
      rz=refc(3)
      call dshap(fbec8g25,refc,shpr,3,8,1)
c .... shape function and their derivatives computation
c .... compute partial derivatives by centered difference
c .... which is in the file ccshap.for of FEPG library
      return
      end

      real*8 function fbec8g25(refc,n)
c .... shape function caculation
      implicit real*8 (a-h,o-z)
      common /ccbec8g2/ xa(8),ya(8),za(8),una(8),
     &vna(8),wna(8)
      common /vbec8g2/ rctr(3,3),crtr(3,3),coefd(3,9),coefc(3,9)
      dimension refc(3)
      common /coord/ coor(3),coora(27,3)
      x=coor(1)
      y=coor(2)
      z=coor(3)
      rx=refc(1)
      ry=refc(2)
      rz=refc(3)
      goto (1,2,3,4,5,6,7,8) n
1     fbec8g25=+(+1.-rx)/2.*(+1.-ry)/2.*(+1.-rz)/2. 
      goto 1000
2     fbec8g25=+(+1.+rx)/2.*(+1.-ry)/2.*(+1.-rz)/2. 
      goto 1000
3     fbec8g25=+(+1.+rx)/2.*(+1.+ry)/2.*(+1.-rz)/2. 
      goto 1000
4     fbec8g25=+(+1.-rx)/2.*(+1.+ry)/2.*(+1.-rz)/2. 
      goto 1000
5     fbec8g25=+(+1.-rx)/2.*(+1.-ry)/2.*(+1.+rz)/2. 
      goto 1000
6     fbec8g25=+(+1.+rx)/2.*(+1.-ry)/2.*(+1.+rz)/2. 
      goto 1000
7     fbec8g25=+(+1.+rx)/2.*(+1.+ry)/2.*(+1.+rz)/2. 
      goto 1000
8     fbec8g25=+(+1.-rx)/2.*(+1.+ry)/2.*(+1.+rz)/2. 
      goto 1000
1000  return
      end

      subroutine bec8g26(refc,shpr,rctr,crtr)
c .... compute shape functions and their partial derivatives
c .... shapr ---- store shape functions and their partial derivatives
      implicit real*8 (a-h,o-z)
      dimension refc(3),shpr(8,4),rctr(3,3),crtr(3,3)
      external fbec8g26
      rx=refc(1)
      ry=refc(2)
      rz=refc(3)
      call dshap(fbec8g26,refc,shpr,3,8,1)
c .... shape function and their derivatives computation
c .... compute partial derivatives by centered difference
c .... which is in the file ccshap.for of FEPG library
      return
      end

      real*8 function fbec8g26(refc,n)
c .... shape function caculation
      implicit real*8 (a-h,o-z)
      common /ccbec8g2/ xa(8),ya(8),za(8),una(8),
     &vna(8),wna(8)
      common /vbec8g2/ rctr(3,3),crtr(3,3),coefd(3,9),coefc(3,9)
      dimension refc(3)
      common /coord/ coor(3),coora(27,3)
      x=coor(1)
      y=coor(2)
      z=coor(3)
      rx=refc(1)
      ry=refc(2)
      rz=refc(3)
      goto (1,2,3,4,5,6,7,8) n
1     fbec8g26=+(+1.-rx)/2.*(+1.-ry)/2.*(+1.-rz)/2. 
      goto 1000
2     fbec8g26=+(+1.+rx)/2.*(+1.-ry)/2.*(+1.-rz)/2. 
      goto 1000
3     fbec8g26=+(+1.+rx)/2.*(+1.+ry)/2.*(+1.-rz)/2. 
      goto 1000
4     fbec8g26=+(+1.-rx)/2.*(+1.+ry)/2.*(+1.-rz)/2. 
      goto 1000
5     fbec8g26=+(+1.-rx)/2.*(+1.-ry)/2.*(+1.+rz)/2. 
      goto 1000
6     fbec8g26=+(+1.+rx)/2.*(+1.-ry)/2.*(+1.+rz)/2. 
      goto 1000
7     fbec8g26=+(+1.+rx)/2.*(+1.+ry)/2.*(+1.+rz)/2. 
      goto 1000
8     fbec8g26=+(+1.-rx)/2.*(+1.+ry)/2.*(+1.+rz)/2. 
      goto 1000
1000  return
      end

      subroutine bec8g27(refc,shpr,rctr,crtr)
c .... compute shape functions and their partial derivatives
c .... shapr ---- store shape functions and their partial derivatives
      implicit real*8 (a-h,o-z)
      dimension refc(3),shpr(8,4),rctr(3,3),crtr(3,3)
      external fbec8g27
      rx=refc(1)
      ry=refc(2)
      rz=refc(3)
      call dshap(fbec8g27,refc,shpr,3,8,1)
c .... shape function and their derivatives computation
c .... compute partial derivatives by centered difference
c .... which is in the file ccshap.for of FEPG library
      return
      end

      real*8 function fbec8g27(refc,n)
c .... shape function caculation
      implicit real*8 (a-h,o-z)
      common /ccbec8g2/ xa(8),ya(8),za(8),una(8),
     &vna(8),wna(8)
      common /vbec8g2/ rctr(3,3),crtr(3,3),coefd(3,9),coefc(3,9)
      dimension refc(3)
      common /coord/ coor(3),coora(27,3)
      x=coor(1)
      y=coor(2)
      z=coor(3)
      rx=refc(1)
      ry=refc(2)
      rz=refc(3)
      goto (1,2,3,4,5,6,7,8) n
1     fbec8g27=+(+1.-rx)/2.*(+1.-ry)/2.*(+1.-rz)/2. 
      goto 1000
2     fbec8g27=+(+1.+rx)/2.*(+1.-ry)/2.*(+1.-rz)/2. 
      goto 1000
3     fbec8g27=+(+1.+rx)/2.*(+1.+ry)/2.*(+1.-rz)/2. 
      goto 1000
4     fbec8g27=+(+1.-rx)/2.*(+1.+ry)/2.*(+1.-rz)/2. 
      goto 1000
5     fbec8g27=+(+1.-rx)/2.*(+1.-ry)/2.*(+1.+rz)/2. 
      goto 1000
6     fbec8g27=+(+1.+rx)/2.*(+1.-ry)/2.*(+1.+rz)/2. 
      goto 1000
7     fbec8g27=+(+1.+rx)/2.*(+1.+ry)/2.*(+1.+rz)/2. 
      goto 1000
8     fbec8g27=+(+1.-rx)/2.*(+1.+ry)/2.*(+1.+rz)/2. 
      goto 1000
1000  return
      end

      subroutine bec8g28(refc,shpr,rctr,crtr)
c .... compute shape functions and their partial derivatives
c .... shapr ---- store shape functions and their partial derivatives
      implicit real*8 (a-h,o-z)
      dimension refc(3),shpr(8,4),rctr(3,3),crtr(3,3)
      external fbec8g28
      rx=refc(1)
      ry=refc(2)
      rz=refc(3)
      call dshap(fbec8g28,refc,shpr,3,8,1)
c .... shape function and their derivatives computation
c .... compute partial derivatives by centered difference
c .... which is in the file ccshap.for of FEPG library
      return
      end

      real*8 function fbec8g28(refc,n)
c .... shape function caculation
      implicit real*8 (a-h,o-z)
      common /ccbec8g2/ xa(8),ya(8),za(8),una(8),
     &vna(8),wna(8)
      common /vbec8g2/ rctr(3,3),crtr(3,3),coefd(3,9),coefc(3,9)
      dimension refc(3)
      common /coord/ coor(3),coora(27,3)
      x=coor(1)
      y=coor(2)
      z=coor(3)
      rx=refc(1)
      ry=refc(2)
      rz=refc(3)
      goto (1,2,3,4,5,6,7,8) n
1     fbec8g28=+(+1.-rx)/2.*(+1.-ry)/2.*(+1.-rz)/2. 
      goto 1000
2     fbec8g28=+(+1.+rx)/2.*(+1.-ry)/2.*(+1.-rz)/2. 
      goto 1000
3     fbec8g28=+(+1.+rx)/2.*(+1.+ry)/2.*(+1.-rz)/2. 
      goto 1000
4     fbec8g28=+(+1.-rx)/2.*(+1.+ry)/2.*(+1.-rz)/2. 
      goto 1000
5     fbec8g28=+(+1.-rx)/2.*(+1.-ry)/2.*(+1.+rz)/2. 
      goto 1000
6     fbec8g28=+(+1.+rx)/2.*(+1.-ry)/2.*(+1.+rz)/2. 
      goto 1000
7     fbec8g28=+(+1.+rx)/2.*(+1.+ry)/2.*(+1.+rz)/2. 
      goto 1000
8     fbec8g28=+(+1.-rx)/2.*(+1.+ry)/2.*(+1.+rz)/2. 
      goto 1000
1000  return
      end

      subroutine tbec8g2(refc,coor,coorr,coefr,rc)
c .... compute coordinate value and Jacobi's matrix rc
c .... by reference coordinate value
      implicit real*8 (a-h,o-z)
      dimension refc(3),coor(3),coorr(3,8),coefr(8,3),rc(3,3)
      common /ccbec8g2/ x(8),y(8),z(8),un(8),vn(8),wn(8)
      external ftbec8g2
      do 100 n=1,8
      x(n)=coorr(1,n)
      y(n)=coorr(2,n)
      z(n)=coorr(3,n)
100   continue
      do 200 n=1,8
      un(n)=coefr(n,1)
      vn(n)=coefr(n,2)
      wn(n)=coefr(n,3)
200   continue
      rx=refc(1)
      ry=refc(2)
      rz=refc(3)
      call dcoor(ftbec8g2,refc,coor,rc,3,3,1)
c .... coordinate value and their partial derivatives caculation
c .... compute partial derivatives by centered difference
c .... which is in the file ccshap.for of FEPG library
      return
      end

      real*8 function ftbec8g2(refc,n)
c .... coordinate transfer function caculation
      implicit real*8 (a-h,o-z)
      dimension refc(3)
      common /ccbec8g2/ x(8),y(8),z(8),un(8),vn(8),wn(8)
      common /vbec8g2/ rctr(3,3),crtr(3,3),coefd(3,9),coefc(3,9)
      rx=refc(1)
      ry=refc(2)
      rz=refc(3)
      goto (1,2,3) n
1     ftbec8g2=+(+(+1.-rx)/2.*(+1.-ry)/2.*(+1.-rz)/2.)*x(1)
     & +(+(+1.+rx)/2.*(+1.-ry)/2.*(+1.-rz)/2.)*x(2)+(+(+1.+
     & rx)/2.*(+1.+ry)/2.*(+1.-rz)/2.)*x(3)+(+(+1.-rx)/2.*(+
     & 1.+ry)/2.*(+1.-rz)/2.)*x(4)+(+(+1.-rx)/2.*(+1.-ry)/2.*
     & (+1.+rz)/2.)*x(5)+(+(+1.+rx)/2.*(+1.-ry)/2.*(+1.+rz)/
     & 2.)*x(6)+(+(+1.+rx)/2.*(+1.+ry)/2.*(+1.+rz)/2.)*x(7)+(+
     & (+1.-rx)/2.*(+1.+ry)/2.*(+1.+rz)/2.)*x(8)
      goto 1000
2     ftbec8g2=+(+(+1.-rx)/2.*(+1.-ry)/2.*(+1.-rz)/2.)*y(1)
     & +(+(+1.+rx)/2.*(+1.-ry)/2.*(+1.-rz)/2.)*y(2)+(+(+1.+
     & rx)/2.*(+1.+ry)/2.*(+1.-rz)/2.)*y(3)+(+(+1.-rx)/2.*(+
     & 1.+ry)/2.*(+1.-rz)/2.)*y(4)+(+(+1.-rx)/2.*(+1.-ry)/2.*
     & (+1.+rz)/2.)*y(5)+(+(+1.+rx)/2.*(+1.-ry)/2.*(+1.+rz)/
     & 2.)*y(6)+(+(+1.+rx)/2.*(+1.+ry)/2.*(+1.+rz)/2.)*y(7)+(+
     & (+1.-rx)/2.*(+1.+ry)/2.*(+1.+rz)/2.)*y(8)
      goto 1000
3     ftbec8g2=+(+(+1.-rx)/2.*(+1.-ry)/2.*(+1.-rz)/2.)*z(1)
     & +(+(+1.+rx)/2.*(+1.-ry)/2.*(+1.-rz)/2.)*z(2)+(+(+1.+
     & rx)/2.*(+1.+ry)/2.*(+1.-rz)/2.)*z(3)+(+(+1.-rx)/2.*(+
     & 1.+ry)/2.*(+1.-rz)/2.)*z(4)+(+(+1.-rx)/2.*(+1.-ry)/2.*
     & (+1.+rz)/2.)*z(5)+(+(+1.+rx)/2.*(+1.-ry)/2.*(+1.+rz)/
     & 2.)*z(6)+(+(+1.+rx)/2.*(+1.+ry)/2.*(+1.+rz)/2.)*z(7)+(+
     & (+1.-rx)/2.*(+1.+ry)/2.*(+1.+rz)/2.)*z(8)
      goto 1000
1000  return
      end

      subroutine ebec8g2(refc,coef,coorr,coefr,coefd)
c .... compute coef value and their partial derivatives
c .... by reference coordinate value
      implicit real*8 (a-h,o-z)
      dimension refc(3),coef(3),coorr(3,8),coefr(8,3),coefd(3,3)
      external febec8g2
      rx=refc(1)
      ry=refc(2)
      rz=refc(3)
      call dcoef(febec8g2,refc,coef,coefd,3,3,2)
c .... coef value and their partial derivatives caculation
c .... compute partial derivatives by centered difference
c .... which is in the file ccshap.for of FEPG library
      return
      end

      real*8 function febec8g2(refc,n)
c .... coef function caculation
      implicit real*8 (a-h,o-z)
      dimension refc(3)
      common /ccbec8g2/ xa(8),ya(8),za(8),un(8),vn(8),wn(8)
      common /vbec8g2/ rctr(3,3),crtr(3,3),coefd(3,9),coefc(3,9)
      common /coord/ coor(3),coora(27,3)
      x=coor(1)
      y=coor(2)
      z=coor(3)
      rx=refc(1)
      ry=refc(2)
      rz=refc(3)
      goto (1,2,3) n
1     febec8g2=+(+(+1.-rx)/2.*(+1.-ry)/2.*(+1.-rz)/2.)*un(1)
     & +(+(+1.+rx)/2.*(+1.-ry)/2.*(+1.-rz)/2.)*un(2)+(+(+1.+
     & rx)/2.*(+1.+ry)/2.*(+1.-rz)/2.)*un(3)+(+(+1.-rx)/2.*
     & (+1.+ry)/2.*(+1.-rz)/2.)*un(4)+(+(+1.-rx)/2.*(+1.-ry)/
     & 2.*(+1.+rz)/2.)*un(5)+(+(+1.+rx)/2.*(+1.-ry)/2.*(+1.+
     & rz)/2.)*un(6)+(+(+1.+rx)/2.*(+1.+ry)/2.*(+1.+rz)/2.)*un(7)
     & +(+(+1.-rx)/2.*(+1.+ry)/2.*(+1.+rz)/2.)*un(8)
      goto 1000
2     febec8g2=+(+(+1.-rx)/2.*(+1.-ry)/2.*(+1.-rz)/2.)*vn(1)
     & +(+(+1.+rx)/2.*(+1.-ry)/2.*(+1.-rz)/2.)*vn(2)+(+(+1.+
     & rx)/2.*(+1.+ry)/2.*(+1.-rz)/2.)*vn(3)+(+(+1.-rx)/2.*
     & (+1.+ry)/2.*(+1.-rz)/2.)*vn(4)+(+(+1.-rx)/2.*(+1.-ry)/
     & 2.*(+1.+rz)/2.)*vn(5)+(+(+1.+rx)/2.*(+1.-ry)/2.*(+1.+
     & rz)/2.)*vn(6)+(+(+1.+rx)/2.*(+1.+ry)/2.*(+1.+rz)/2.)*vn(7)
     & +(+(+1.-rx)/2.*(+1.+ry)/2.*(+1.+rz)/2.)*vn(8)
      goto 1000
3     febec8g2=+(+(+1.-rx)/2.*(+1.-ry)/2.*(+1.-rz)/2.)*wn(1)
     & +(+(+1.+rx)/2.*(+1.-ry)/2.*(+1.-rz)/2.)*wn(2)+(+(+1.+
     & rx)/2.*(+1.+ry)/2.*(+1.-rz)/2.)*wn(3)+(+(+1.-rx)/2.*
     & (+1.+ry)/2.*(+1.-rz)/2.)*wn(4)+(+(+1.-rx)/2.*(+1.-ry)/
     & 2.*(+1.+rz)/2.)*wn(5)+(+(+1.+rx)/2.*(+1.-ry)/2.*(+1.+
     & rz)/2.)*wn(6)+(+(+1.+rx)/2.*(+1.+ry)/2.*(+1.+rz)/2.)*wn(7)
     & +(+(+1.-rx)/2.*(+1.+ry)/2.*(+1.+rz)/2.)*wn(8)
      goto 1000
1000  return
      end

